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In massive gravity, galileon, and braneworld explanations of cosmic acceleration, force modifica- 
tions are screened by nonlinear derivative self-interactions of the scalar field mediating that force. 
Interactions between the field of a central body ( "A" ) and an orbiting body ( "_B" ) imply that body 
B does not move as a test body in the field of body A if the orbit is smaller than the Vainshtein 
radius of body B. We find through numerical solutions of the joint field at the position of B that 
the A-field Laplacian is nearly perfectly screened by the B self-field, whereas first derivative or net 
forces are reduced in a manner that scales with the mass ratio of the bodies as {Mb / MaY^^ ■ The 
latter causes mass-dependent reductions in the universal perihelion precession rate, with deviations 
for the Earth-Moon system at the ~ 4% level. In spite of universal coupling, which preserves the 
microscopic equivalence principle, the motion of macroscopic screened bodies depends on their mass 
providing in principle a means for testing the Vainshtein mechanism. 

PACS numbers: 98.80.-k 



I. INTRODUCTION 



The current acceleration of the cosmic expansion is one of the most puzzling aspects of modern cosmology. Aside 
from a cosmological constant whose smallness remains unexplained, the simplest models typically involve an additional 
scalar field either implicitly or explicitly. Universal coupling of this field to matter would produce gravitational strength 
fifth forces and naively be excluded by laboratory and solar system bounds. Viable models must therefore implement 
a so-called screening mechanism to hide fifth forces locally. 

Screening mechanisms invoke nonlinearity in the field equations to change the nature of the fifth force in high density 
regions. For example, the chameleon mechanism increases the mass of the field in deep gravitational potentials [1] 
whereas the symmetron mechanism changes its coupling to matter [5] . A third possibility is the Vainshtein mechanism 
[31 0] , first introduced in the context of massive gravity to suppress the propagation of additional helicity modes [51 [5] . 
Here nonlinear derivative interactions of the field act to screen the fifth force within the so-called Vainshtein radius 
around a matter source. The Vainshtein mechanism occurs not only in modern incarnations of Boulware-Deser [7] 
ghost-free massive gravity [SHU] but also in galileon cosmology [T^HTB] and braneworld models. Indeed it is in the 
braneworld model of Dvali, Gabadadze and Porrati (DGP) that it has been best studied [2D1EB] . 

Interestingly, these mechanisms are distinguished by how screened bodies fall in external fields |27j . As a consequence 
of universal coupling, all unscreened test bodies fall in the same way and obey a microscopic equivalence principle. In 
the chameleon and symmetron models, screened bodies do not respond to external fields. In the Vainshtein mechanism 
they do, but only if those fields have wavelengths long compared to the Vainshtein radius. These differences arise 
because the self-field of the screened body and the external field do not in general superimpose but rather interfere 
in a manner dependent on the nonlinear interaction |28) . 

In this paper, we consider the Vainshtein mechanism in the near-field limit. In particular we study the motion of 
two bodies which are separated by much less than their individual Vainshtein radii and look for apparent violations 
of the equivalence principle. The two-body problem is particularly relevant since it has been shown that for the orbit 
of test bodies, there is a universal anomalous precession rate induced by a Vainshtein-screened scalar field that is 
potentially measurable in next generation solar system tests of General Relativity [2T1 [5S] . However in the Earth- 
Moon system the Moon is screened on scales out to nearly a parsec for cosmologically motivated models and cannot 
be considered as a test body in the Earth's field. In principle, this nonlinearity can affect the interpretation of lunar 
ranging tests for anomalous precession and more generally lead to results that depend on the nature of the orbiting 
body. To understand this system, the field must be solved jointly in the presence of both sources. For definiteness, we 
will implement the Vainshtein mechanism in the DGP model but our results readily apply to the galileon and massive 
gravity incarnations as well. 

The remainder of the paper is organized as follows. We briefly review the Vainshtein mechanism and the spherically- 
symmetric one-body solution in Sec.|ll] In Sec. |III[ we discuss violation of the superposition principle in the two-body 



case. We present numerical results and their scaling with the two-body parameters in Sec. IV Details of the numerical 
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scheme is given in Appendix \K\ We discuss the imphcations of these results in Sec. [Vj 

II. VAINSHTEIN MECHANISM 
A. DGP example 

As an example of models that accommodate the Vainshtein mechanism, we consider the DGP braneworld model 
[TO] . In the DGP model, there is a dynamical degree of freedom representing the bending of the brane embedded in 
the five-dimensional bulk which we denote '/>(x, t). In the quasi-static limit where its time-derivatives can be neglected, 
its equation of motion becomes [2TJ [23] 

3/3(t)^0+ij7V[0,(/.] -8^G5p, (1) 

where 

m - 1 ± 2Hr, (^1 + , (2) 

N[<t)A, ^b] - rliV^cpAV^^B ~ V,V,<t>AW'^B), (3) 

and 5p = p — p, the density fluctuation from the cosmic mean. Spatial derivatives here are in comoving coordinates. 
The key parameter in this model is the crossover scale 



2G" 



(4) 



the ratio between the five dimensional Newton constant G5 and the four dimensional one G. 

The -I- sign in I3{t) corresponds to the normal branch solution while the — sign corresponds to the self-accelerating 
solution. In this paper, we only consider the solutions with /3 > since /3 < is associated with ghost instabilities 
|30j . Furthermore, since we are interested in static solutions at the current epoch a = 1, we set /3 = 1. All results can 
be rescaled to the general /3 case by replacing Tc with Vc/y/P and densities 5p with 5p/ j3, or equivalently the masses 
of all bodies. 

We have written the nonlinear operator N[(j)A, 0b] in bilinear form allowing for two separate fields since in the two 
body calculation that follows it will be useful to consider the interference between the two individual fields. It is the 
nonlinearity of this operator that is responsible for both the Vainshtein mechanism and the lack of a superposition 
principle for the brane bending mode. 

B. One Body Solution 

The Vainshtein mechanism and the scale associated with it can be illustrated with simple analytic one-body so- 
lutions. For a spherically symmetric object with a top-hat constant density, the scalar field equation ([T]) reduces 
to 
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where the top-hat density profile is 



2 / d<j)V A (f(j) dcf) 



i?2 \dR RdR^ dR 



= SnGSp, (5) 



R> 



with Tg as the radius of the source. Eq. ([5]) can be integrated by multiplying both sides by R^dR , resulting in 

3R^^ + 2Rrl(^X^8nG / ^i?' ^'^''^f/' " (7) 
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As this is a quadratic equation for d(f>/dR we immediately obtain 

d(j) _ 3R jg{r,) R < r„ 



X 



dR 4r2 ]g{R) R> 
where 



(8) 



2\ 1/3 



and the Schwarzschild radius 

= 2GM = ^Sporl (10) 



The radius is caUed the Vainshtein radius. For Ts < i? ^ r*, we obtain g{R) oc so c>c i?^/^ +const. 

This means that for test bodies the correction to Newtonian forces around the source vanishes in the hmit ^ r*, 
and thus Einstein gravity is recovered (35]. This condition is satisfied if Vg ^ 9r^/(8r^) which is the case for stars and 
even the Earth and the Moon as fong as we set Tc to be a cosmological scale. Equivalently, for spherically symmetric 
systems the Vainshtein mechanism is active whenever the mean enclosed overdensity 3M/47ri?"^ is greater than 
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where pcr is the critical density. 

We can in fact obtain the full solution for in closed form For i? > r^, defining x = R/r^, and A = 3rJ/(4r^), 
we can express the solution of Eq. ([s]) in terms of the hypergeometric function, 

0ox(a;) = Ci + -^fi^), for x > x^, (12) 



where Ci is an integral constant, Xs = rs/r^, and 



fix) EE x' 



(13) 



On the other hand, in the internal region, R < r^, the solution of Eq. (|8|) is 



</'in(a;) = —x^g{xs) + C2, for x < Xs- (14) 



Note that in this case, the two pieces of the nonlinear terms combine and imply 



More generally, for any (j)A 



[(V2(/)i„)2 „ (V,V,0in)'] = ^(V20in)2. (15) 



N[^A, 0in] = ^rlV^C^AV^^in- (16) 



This relation is specific to the interior of a top hat (f)ia but nonetheless will be useful when approximating the nonlinear 
term. 

The inner and external solutions should be continuous at x = a;^ . Without loss of generality, we can take Ci — 
due to the shift symmetry of the scalar field equation of motion. Then we obtain the solution for a single source, 

4>{x) = (l)in{x)Q{Xs — x) + (/)cxix)Q{x — Xs) with 

^M^^fix), cbi^ix) = ^ [f{xs) + gixs)ix^ - xD] . (17) 
For R and r^, these solutions can be approximated as 

^ ^ „2 , Sry^rl/' r[-2/3]r[7/6] 

0in i?) ~ —p — jT^R + ^77i? — + ^ T ' (18) 
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with a Laplacian of 



0^1/2 o 1/2 



That the radial dependence of these solutions is proportional to oc M^/^ is an indication that mass sources do 
not linearly superimpose within the Vainshtein radius. On the other hand in the opposite limit of i? 3> 3> , 



3i? V" 18^0 

Leading order linearity in cx M implies that mass sources do superimpose in this limit. Since the leading order 
term in the Laplacian — 0, residual effects go as 

VV(i?)-^i^- (21) 

These approximate forms will be useful in constructing scaling relations and boundary conditions for the two body 
problem. 

III. TWO BODY PROBLEM 

The brane bending field of a single body derived in the previous section suffices to study the motion of test particles 
around that body. For test particles, there is a universal anomalous precession rate, dependent only on the crossover 
scale Tc, that can be used to test the Vainshtein mechanism [5TJ However for most realistic orbiting bodies, 
their own Vainshtein radius is too large for them to be considered test particles. This is in particular true for the 
Earth-Moon system where the orbit of the Moon is well within its own Vainshtein radius. 

Through the nonlinearity of the Vainshtein effect, the orbiting body's self- field interferes with that of the central 
body and can in principle affect its orbit. Indeed the scaling of the single body field as \fM is an indication of the 
non-superimposability of solutions within the Vainshtein radius. Two sources will add as \JMa + Mb for distances 
from the center of mass much larger than the separation. 

In this section, we consider a two-body problem such as the Earth-Moon system to study non-superimposability 



of solutions. We begin in ^III A with the parameterization of the two body system in terms of the physical scales in 
the problem. We describe how screening operates directly on second derivatives of the field and indirectly on first 
derivatives, or average forces in §111 B[ We examine the geometry of screening in |III C| and introduce our screening 



statistics and their scaling properties in \ III D 



A. Model Parameters 

Given spherically symmetric masses, the system has axial symmetry and so we use cylindrical coordinates (r, 0, z). 
We assume that the two bodies, denoted as A and i3, are separated from each other by c?, and they are located at 
(r, z) = (0, 0) for body A and (r, z) = (0, —d) for body B. We denote their respective Schwarzschild and physical radii 
as rgA,rgB and rsAi^sB- In the following we assume body A is heavier than body B. The schematic illustration of 
the two-body set up is shown in Fig. [l] 

In our numerical solutions below, we choose parameters that reflect the Earth ( "E" ) and Moon ( "M" ) where possible. 
In the actual Earth-Moon system d = 3.8 x 10^ km and choosing ~ cHq^ as appropriate for cosmologically- motivated 
models would require an 10^ dynamic range between the Vainshtein radius of the Earth and the separation d since 



Furthermore between the separation scale d and the body radii are several more orders of magnitude 

rsE = 6378 km = 1.7 x lO'^d, 

rsM = 1737 km = 4.6 x lO'^d. (23) 
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FIG. 1: Schematic illustration of the two body problem (not drawn to scale). Physical radii are shown by the solid lines whereas 
the Vainshtein radii, r,A and r*s are shown in dashed lines. Boundary conditions are set at r = L and z — ±L with vanishing 
deviations from superposition of A and B. 



On the other hand, we shall see that the most important properties of the two body solution can be expressed as a 
fmiction of the mass ratio Mb /Ma, whereas other properties can be inferred by examining the scahngs of results with 
rc/d and rg/d. 

Our fiducial choice will therefore be to take d and masses, or equivalently Schwarzschild radii, from the Earth-Moon 
system 



rgA _ ^GMe 
d ~ d 

rgB ^ ^ GMm 
d ~ d 



2.3 X 10"", 

rgA f Mm \ 
d [MeJ 



(24) 



Since we are interested in effects around body B we fix TsA/d = 0.3 and explore the dependence on rss- Likewise we 
explore the dependence of results on rdd. Our fiducial choices for these parameters are listed in Table. |lj 



Mb /Ma 


1/80 


rgA/d 


2.3 X 10"" 


Tc/d 


10* 


TaA/d 


0.3 


Tss/d 


0.1 


TgB/d 


2.8 X 10~" 


r*A/d 


58.91 


r^s/d 


13.67 


L/d 


286.7 


h/d 


0.0125 



TABLE 1; The fiducial parameters for the two-body problem where Mb /Ma and VgA/d mimic the Earth-Moon system. The 
top set represents the primary parameters whose individual variation we explore below. The middle set represents derived 
parameters. The last two parameters are numerical and give the box and grid size (see Appendix M. 



It is useful both numerically and analytically to express the two body problem in terms of deviations from the 
superposition principle. Denoting the full solution as 

= 0A + 0B + 0A, (25) 

where (\)a and 0s are the solutions to the single body equations, the field equation for the deviation is 

3V20A + iV[0A, (/-a] + 2iV[0^ + 0B, 0a] = -27V[0A, 0b]. (26) 
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The interference term N[(I)a, 0b] can be considered as the source of deviations from superposition. 
We solve Eq. ( 26 ) numerically with boundary conditions 

<Pa{L, z) = 0A(r, -i) = 0a(?', L) = 0, 



(27) 



corresponding to the case L ^ r*^ where we can superpose single-body solutions. Finally the equations are solved 
by finite difference on an inhomogeneous grid with minimum spacing h/d (see Appendix [A|). 



B. Screening 



Before turning to numerical solutions of the two body system it is useful to examine the structure of Eq. ( 26 ) and 



develop an analytic intuition for the results. In particular, Eq. (26) admits screening solutions where the field of body 



B suppresses some aspect of the field of body A around itself and vice versa. 

To see this screening behavior note that both (j)^ = —4>a and 0a = —ipB are solutions to this equation outside of 
the sources where 3V^(/>^.b + N[(I)a,b, 'Pa,b] = 0. In particular, around body B we might expect (/)a = — such that 
it cancels the effect of body A. However these source-free solutions would not match the boundary conditions at the 
sources themselves and so what occurs in reality is more complicated. 

Screening really occurs in the second derivatives of the field. Note that we are also free to add a pure gradient to 
any of these source-free solutions, e.g. 



D 



(28) 



and still solve the equations. Thus we might expect that screening operates by replacing (j)A with a version of itself 
that is linearized over some region of influence around body B that can be much larger than the physical size of the 
body (cf. [22). 



We start with the simple expectations from approximating the nonlinear term with the Laplacian in Eq. (16). In 
this approximation the more general screening expectation becomes 



V^^A 



V^(j)B > V20^ 
V20B < V20^ 



(29) 



with the constants C and D providing the appropriate matching of the regimes. Note that if (j)A is already nearly 
linear around body i3, as is the case for a distant source [27] we expect no self- field effect on the motion of body B. 
On the other hand we know that there must be a near field effect on forces between the bodies: without screening 
the force from body A on body B 



Fab 'X Mb V Ma, 



(30) 



whereas that from body B on body A 



Fba oc -MaVMb- 



(31) 



This violation of Newton's third law would cause momentum non-conservation in the joint system. 

To get a rough sense for the size and scaling of near field effects we can replace cjyA oc y/R with a linearization of 
itself across the region on the z-axis where V^(j)B > V^0a 



1 

l + (AfB/M^)i/3' 



•2inax 



1-(Mb/M4)1/3' 



(32) 



0a(O, Z) -0^(0, Z) H (Z - Zmin) + 0a(O, Zmin)- (33) 

Zm ax ZjTi 1 n 



In this crude approximation, the gradient 920^(0, —rf) is reduced by 

2/3 



dz(l)A ,r, j\ 3 / Mb 



[Mb^"' 



O 

\Ma 



(34) 



0.0 



0.5 



1.0 



1.5 



2.0 



r/d 



FIG. 2: Two-body nonlinearity source function F[r/d,z/d\ oc N[(I)a,4'b\ with A at (0,0) and B at (0,-1) with parameters of 
the fiducial model (see Tab.|l|. The region interior to the bodies, within the semicircles, is not shown. 

independently of the physical size of the bodies and the value of Tc/d ^ 1. For the Earth-Moon mass ratio this is a 
^ 2% correction of the gradient even though the second derivative is screened across a much larger range, Az/d ^ 0.5, 
than the physical size of the Moon Az/d ^ 0.005. Nonetheless the dependence on the mass ratio represents an apparent 
violation of the equivalence principle. 

Note that momentum conservation in the joint system would imply that at body A there is near complete screening 
of the force from body B 



to determine how well Newton's third law Fba/Fab = — 1 is satisfied. 

Note that in the opposite limit d <C r <C rA* , where the sources add as V Ma + Mb , we know that the field of body 
A screens that of body B more directly. The (p^ field here is simply the difference between the joint and individual 
sources 




(35) 



for TsA, ^sB d. As a check of our numerical results we will examine 



Fab _ Mb {dz(pA + <9^0a) |o 



(36) 



Fba Ma (92</>B + (9^0a) |o,o 




(37) 



for Mb <C Ma implying 



0A 



4>B- 



(38) 



While these considerations provide a qualitative guide to results, the specific form of the second derivatives in 
N[(j)A,4'B] lead to important consequences for the geometry of the screening around body B which we shall now 
consider. 



C. Toroidal Geometry 



In order to understand the geometry of screening around body B, let us examine the source to the 0a held 



Eq. (26) 



2N[(I)A, (f>B] = -^Tc FVId, z/d], (39) 

where the geometry is determined by a universal function in cylindrical coordinates scaled to the separation d 

F[r,A= (^2+^2)7/4(^2 + +1)2)7/4 ' (^0) 

This function is plotted in Fig. [2j Instead of the roughly spherical geometry that the arguments based on the Laplacian 
would predict, the true interference of the self-field of B on that of A is toroidal around body B. The basic reason 
for this geometry is that along the z axis the cross terms cancel given the difference structure in Eq. ([3|. 
On the other hand, interior to body B, the nonlinear source becomes 



2iV0A,0B = — 773 ] — (41) 

rf/'(r2 + z2)3/4 

and does not vanish for r — Q but rather approaches a constant for d ^ Tsb- Note that for a small physical size 
TsB, r'^ + z'^ w d^ in the interior. The interference term is nearly constant and approximates an effective density 
of y/PBPA,cS where pA.cS = Myi/(47rc?^/3). Unlike the true density, 2N[(I)a, (I)b] does not vanish in the exterior but 
has a jump in value at the r^B boundary. We shall see that this jump causes a similar discontinuity in the second 
derivatives of (pA at the boundary. 

Finally, although 2N[4>a,4>b] peaks at the bodies themselves (saturated in Fig. [2]), this does not necessarily mean 
that the relative impact on the joint field peaks there. The single body Laplacians also peak there and we must 
examine deviations with respect to those fields in the results that follow. 



D. Screening Statistics 



It is therefore useful to introduce our primary, or Laplacian, screening statistic for the relative impact of the field 
of body B on that of A 



Q2{r, z) 



(42) 



If screening is complete, Q2 — >■ —1. 

It is again useful to derive rough analytic scalings for the screening statistic in various limits. When either A 
or B dominate in the second derivatives we expect that the second derivatives of the correction field is small 



compared with the dominant one and screens the subdominant one. We therefore expect the field equation (26 1 to be 
satisfied approximately by dropping terms nonlinear in 0a- Furthermore, the SV^^a terms is small compared with 
the nonlinear terms deep within the Vainshtein radius and so 



N[(t)A + (l'B,(l'A\ 



'N[<l)A,(t>l 



We can further approximate the left hand side of Eq. (|43|) using Eq. (|16) 

N[(f>A + <pB,<l>A] 



)V^ 



OA- 



(43) 



(44) 



This relation is exact only for the interior field and for V'^(t)B ^ '^'^4>A but gives a reasonable heuristic description 
elsewhere. This approximation is useful in that it allows us to solve directly for V^^a, and hence 
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0.5 1 1.5 0.5 1 1.5 

r/d rid 

FIG. 3: Laplacian-screening statistic Qi — V^0a/V^(^a (left: analytic; right: numerical). The red areas indicate a large 
relative deviation from the superposition solution in the Laplacian of the field. In the interior of body B screening of the 
Laplacian of body A is nearly complete whereas in the exterior there is a toroidal region of mutual screening. The analytic 
description is in good quahtative agreement near body B. Parameters are the fiducial choice of Tab. |l] 



The left hand panel in Fig. |3] shows the results of Q2ir, z) obtained from Eqs. (451 and (46 1. 

In fact for small radius r <C rf, we can approximate + ~ d^. Then on the z = —d axis, Q2 is approximately 
given by 



1 [Ma 

2 \Mb 

1+r^) 

\MbJ 



1/2 



1/2 



r \ 3/2 

d 



(47) 
(48) 



Note that saturation to Q2 = — 1 in the interior increases with increasing Mb /Ma and decreasing Tsb as expected 
from the fact that the interior value of V'^(j)B scales with these parameters. We can infer from this scaling that for 
realistic situations where rgB ^ d (see Eq. 23 1, that Q2 ^ —1 for any Mb < Ma- At the body surface, there is a 
jump to a maximum value of —3/2. We shall sec from numerical results that this maximum value is only approximate 
since Eq. (44) is not strictly valid in this limit. 



Nonetheless, the qualitative aspects of Q2 indicate that the impact on forces is as a redistribution of force changes 
across a toroidal region around body B rather than a linearization across a quasi-spherical one. Since the volume 
of the regions are comparable, we expect the scaling behavior of Eq. ( 34 1 to be roughly satisfied. To quantify this 
expectation, we define the force-screening statistic 



Qi{r) 



dz4>/ 



z— — d 



(49) 



Note that along z = —d the gradient of the q 
to Qi. Again if screening of forces from (pA 



'B field is along the r-direction and hence its screening does not contribute 
is complete then Qi —1. 



For distances R = \/r'^ + z'^ from the bodies that are large compared with the separation d, Eq. (371 and (21 1 imply 



Q2 



r, MB 



d-^R-^r, 
R > r^A 



(50) 



which will also be useful in checking our numerical results. In particular, Q2 does not depend explicitly on Tc aside 
from setting the transition scale r^,A■ In the R ^ r*^ regime this independence is due to the vanishing of V^0^ to 
leading order. Other statistics do not share this independence. To see the more general dependence on Tc note that 
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the field equation for deviations from superposition Eq. ( 26 ) only has no explicit Tc dependence when the nonlinear 
terms dominate 

N[(j)A,M + 2N[cj)A + (bs, M ~ -2iV[0^, 0s] (51) 

and hence fractional effects of relative to (j)A or (/jb have no Tc dependence. To see when this approximation is 
valid, take the opposite — limit. In that case, 

« ~2N[(bA, <Pb] oc rl r, ^ 0. (52) 



The linear and non-linear terms in Eq. ( 26 ) become comparable when 

max(r2v,Vj0s,r^V,Vj<?;.A) « 1 (53) 



which occurs for the typical distance d around the bodies when r*^ « d. Once 3> d, Eq. (51) becomes valid and 
all statistics lose their dependence on Vc around the bodies. 

IV. RESULTS 

In this section we present numerical solutions to the full nonlinear equation and scaling relations based on them. 
The details of numerical calculations are presented in Appendix [Xj 

We begin by considering the fiducial parameter choices from Tab. |l] In Fig. [s] (right) we show a 2D contour plot 



of the Laplacian-screening statistic Q2 of Eq. (42). As expected from the analytic argument of the previous section, 
Q2 ~ —1 interior to body B indicating nearly full screening. At the surface of the body, it experiences a jump to 
Q2 < which then extends to a toroidal region around the body. In this toroidal region, where the individual body 



0.005 




0.2 p — ^ ■ ' — ■ ' — ■ ■ — T ■ ■ — T ' 1 0.2 




FIG. 4: Force {Qi, top) and Laplacian {Q2, bottom) screening statistics along z — —d as a function of the crossover scale 
Tc (left) and absolute mass scale r^A = 2GMa (right). Other parameters including Mb /Ma are set to the fiducial choices of 
Tab. [ijhere and in the following figures. Arrows indicate the Vainshtein scale of body B, r^B- Near the body the screening 
statistics are independent of Vc and mass scale if r*s 2> d. 
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0.005 




FIG. 5: Screening statistics as a function of the size of body B Tsb (left) and mass ratio Mb /Ma (right) as in Fig. |4] Oi and 
Q2 are independent of Tsb in the exterior r > rsB and converge to constant values for r < Vsb <C d. Likewise for Mb /Ma, Q2 
behaves according to analytic expectations for r ^ d and both Qi and Q2 scale with the mass ratio in the interior r < r^B ~ 0.1. 



Laplacians are comparable, Q2 quantifies both the screening effect of i? on A and A on B. Comparison with Fig. [3] 
(left) shows that our crude analytic approximation of Eq. ( 48 ) captures many of the qualitative effects around body 
B. 

To isolate the effect of screening of B on A, it is useful to examine the force-screening statistic Qi of Eq. (49). The 
field of B is purely radial around B and Qi measures the change in the gradient along the z-direction. As shown in 
Fig. [4] (top), Qi increases toward body B and then smoothly approaches a constant in the interior of the body. For 
this Earth- Moon like system Qi{0) ~ —0.04 or approximately double the crude expectation from Eq. (34). The more 
dramatic changes in Q2 (see Fig. |4j bottom) reflect changes in the radial structure of the joint field. 

The fiducial parameter choices refiect a crossover scale rc/d that is much smaller than the actual Earth- Moon 
system, a coupling strength (3 set to unity by rescaling the masses rg/d, and body sizes r^/d that are large compared 
with the separation. We therefore next test the dependence of our results on these parameters before turning to the 
central dependence on the mass ratio Mb/Ma- 

In Fig. [4] (left), we show the dependence of the force and Laplacian screening statistics on Vc/d. For ^ 10''d, 
the screening statistics scales strongly with Tc as expected from Eq. (52). This behavior saturates once r^s ^ d. In 
Fig.|4j we show the r^.B value corresponding to Vc with an arrow. Indeed so long as Tc 3> lO^d or equivalently r^,B S> d 
the results for both are insensitive to Vc near body B. As expected from Eq. (50), for large Vc/d there is an interval 
d <C r <C r*yi where Q2 ~ —\/Mb/Ma- For r > r^^, Q2 ~ 2Mb /Ma veri fying that we have taken a sufficiently large 
L that the boundary condition (j)^ = is appropriate (see Appendix A 5 for an explicit test) . 

In Fig. |4] (right), we show the dependence on the absolute mass scale or VgA/d at fixed Mb /Ma and other fiducial 
parameters. Recall that changing the mass scale can also be interpreted as changing the parameter /3 in the original 
field equation ([T]). The only dependence of results on the mass scale is through its effect on the Vainshtein scales and 
in the external far field limit relative to the separation d. 

Next we consider the impact of the physical size of body i?, r^B in Fi g, [s] (left). As Tsb shrinks, Q2 interior to the 
body approaches —1. This is expected from our analytic expression, Eq. (48 ), due to the fact that the maximum value 
that V^(/)B attains is controlled by rsB- Once V^^s ^ V^^^i in the interior we expect results to become independent 
of rsB- For the fiducial parameters, Q2 has not quite saturated whereas in Qi it has almost reached its asymptotic 
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value by rss/d = 0.1. In the exterior of body B, Q2 drops increasingly below —1 as TsbIA decreases. This large 
overshoot is not reflected in Q\ and hence represents the screening of the radial body B field itself. Again since the 
maximum V^^s increases as t^b decreases, V^(/)a increases relative to the constant V'^^a as well for the screening 
of the B field. Since statistics in the exterior of body B do not depend on its size, there is likewise no dependence on 
the size of body A, exterior to A. 

We conclude that in the limit of rdd ^ 1 and rss/d <C 1, the main dependence of the screening variables on the 
system parameters is through the mass ratio Mb /Ma- This dependence is shown in Fig. [5] (right). For Q2 interior 
to body B, raising the ratio increases the efficacy of screening body A as expected from Eq. ( |45| . In the exterior, 
it increases the overshooting effect, also as expected. For Qi, screening in the interior scales strongly with the mass 
ratio as expected from E q. (|34| . In both cases, in the far field limit r ^ d, results scale according to the analytic 
expectations of Eq. (37), (|21[), and Eq. (50) to excellent approximation. 



superposition 
3rd law 



6; 

5-0.6 




0.1 0.15 



FIG. 6: Test of Newton's third law or momentum conservation. If the fields of A and B superimposed, then the ratio of forces 
would be — ^ Mb /Ma (red dashed line) and violate the third law (solid blue line). Numerical results (+ points) show that the 
joint solution restores the third law through violation of the superposition principle 0a 7^ 0. 



As shown in Fig. [6j these results for Qi near body B are exactly what is required to restore Newton's third law and 
momentum conservation [see Eq. (36)]. Had the two fields superimposed, then the force ratio would be Mb /Ma 
whereas our numerical results are consistent with —1. 

While our analytic approximations provide a qualitative description of our results, near body B the geometry of 
the screening inhibits their accuracy. It is therefore useful to quantify Qi and Q2 with empirical fits at the center of 
body B. For Qii the numerical results can be described by a scaling relation similar to Eq. (34) but with a finite size 
correction (see Fig. [7| 



Qi(o) 



-0.56 



Mb\ 
Ma) 



0.6 



1 - 0.13 



( Ma 
I Mb 



1/2 



(^) 



3/2 



(54) 



for Mb /Ma ^0.1 and when the correction in brackets is small. In particular, we can extrapolate this fit to the 
physically interesting limit where Tsb/A ^ 1, <3i(0) ~ — 0.56(AfB/Af^)''-^. 

For Qii our results are consistent with approaching —1 in the physically relevant limit 



(32(0, -d) w -1 + 0.23 



TsB 

d 



0-74 /Mb\ 

[maJ 



-0.32 



(55) 



In practice, the approach is somewhat slower than that predicted by Eq. (45) due primarily to inaccuracy in the 
replacement of second derivatives with Laplacians in Eq. (44). Nonetheless if we extrapolate Eq. (55) to the true 
Earth- Moon system where r^s/d = 0.0046 the correction from full screening Q2 = remains only at the percent 
level. 



V. DISCUSSION 



In this paper we studied the apparent violation of the equivalence principle in a model with the Vainshtein mecha- 
nism where nonlinear derivative interactions of the field act to screen the fifth force. We considered the motion of two 
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FIG. 7: Screening statistics at the center of body B, Qi{0) and Q2(0,— d) for various Tsb and Mb /Ma with other parameters 
held fixed to fiducial values. We show numerical results as points (x) as well as the fitting functions in Eqs. (|54[| and ||55[). 



bodies which are separated by much less than their individual Vainshtein radii such as the Earth-Moon system. In 
such a system, the small body B cannot be considered as a test body in the large body A's field 4>a- The nonlinear 
equation for the field allows for screening solutions where the Laplacian of (pA is screened within a region much larger 
in extent than the physical size of body B. The primary effect is that the Laplacian W'^4>a is fully screened. A crude 
estimate of the residual impact on the first derivatives or forces would suggest a suppression of order (Mb/Ma)^^^ 
where Mb/Ma is the mass ratio of the two bodies. Moreover, in the limit where the Vainshtein radii r^A,B are much 
larger than the separation d between the two bodies, we expect relative screening effects to be independent of the 
crossover scale r^- 

We confirmed these qualitative expectations by solving the joint two body system immerically. Cast in terms of the 
deviation of the true field from the superposition solution 0a = 4'~ 4' A ~ 4'b i in the interior of body B there is nearly 
full screening V^^^a ~ in the limit that the size of body B is much smaller than the separation between the 

bodies. On the other hand, the screening of forces from the large body on the small body depend almost exclusively 
on the mass ratio Mb/Ma- From numerical solutions, we found it is given by Qi(0) ~ — 0.56(AfB/M^)'^-^. 

The fifth force introduces an additional contribution to the anomalous perihelion precession rate. In the DGP 
normal branch, the precession rate was obtained by ignoring the non-superimposability and treating the small body 
(e.g., Moon) as a test body. The precession rate (the angle of perihelion advanve A(pDGP during one orbital period 
P) is universal under this assumption and given by 



A(pDGP 3 1 



p 



8rr 



= 7.91 



Horc 



/iarcsec/yr. 



(56) 



This result needs to be revisited in the light of our finding. As screening operates by replacing the large body's field 
with a linearized version of itself, one might think that this affects the precession rate qualitatively by changing the 
radial dependence of the force law. However, the fractional effect is independent of the separation d and the impact 
of non-superimposability comes from the field at the position of the small body B, not at a fixed r. Then the only 
effect of screening is to reduce the large body's force by the factor of 1 -|- Qi and we have a proportional change in 
the precession, 

P " 8r, 



-(i + Qi). 



(57) 
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Precession thus depends on the mass ratio of the bodies and is not universal Different mass bodies will precess at 
different rates leading to an apparent equivalence principle violation. On the other hand, the mass ratio scaling of 
the equivalence principle violation implies that for typical systems the effect will be small. For example, for the mass 
ratio of the Earth-Moon system Mb /Ma — 1/80, leading to a small deviation (4%) from the universal precession 
rate. Nonetheless, in principle the Vainshtein mechanism can be tested by precision tests of the perihelion precession 
of different mass objects. 
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Appendix A: Numerical Techniques 
1. Cylindrical Coordinates 



The field equation for the two-body deviation from superposition cj)/^ given by Eq. ( 26 1 has cylindrical symmetry 
along the 2;-axis separating the bodies (see Fig. [T]) and so 



2r| 
3 



1 



92 



'JAB 



2r| 
3 



AB 



92 



5A 



AB 



Q^2 



dz 



Qj.2 

r dr 



2rl 



1 d(t>AB 
r dr 

92 0A 



92 



^AB 



Qj.2 



A 



'AB 



3 drdz drdz 



where 



and 



>AB 



4>B and we have used 



rf. r 



d(j)a d'^(j)b 
dr dr^ 

d^^a ^ 1 
Qj.2 J, 



Qj.2 Qj. 

d^6b 



dr J dz^ 



^ d'<Pa d^b 

drdz drdz 

d^ rd^b ^ 

dz"^ \ dr"^ 



Id^ 
r dr 



-N[(pA,(l)i 



(Al) 



(A2) 



(A3) 



where (l)a,4'b G 0b, 0ab, 0a}- We treat rgA/d, Mb/Ma^ rsA/d, ^ssjd and r^jd as the model parameters as 
mentioned in Sec. |III A| and so our fundamental unit of length is the separation d. 

In our numerical computations, we impose the regularity condition at r = 0, given as d<j)/dr = 0. As for the other 
boundaries, we assume that superposition holds; namely, we set 0a = at the exterior boundaries. We study the 
effects of the finite computational domain on the numerical solution in Appendix |A 5[ 



2. Nonlinear Coordinates 



To justify the assumption, 0a = at the boundaries, we have to use a box size L 3> r^^A- However, we also want 

3> d and so it is quite difficult to simulate these conditions in the (r, z) coordinate choice. 
Since high resolution is only required near the bodies, we can extend the dynamic range with the help of a nonlinear 
transformation of the radial and axial variables 



Ip\ ^ = c+^c^ (A4) 
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where a is a constant controlling the growth of the grid spacing at large distance. According to these equations, the 
spatial intervals in the original coordinate are given by 

^ ^ A-\p) = {l + ap'), 

^ ^ 5-1(0 = (l + aC^), (A5) 

becoming large at large distance for fixed Ap and relative to the finer resolution at the origin. If we take a = 0.1, 
p = 20d corresponds to r ss 286. TcZ. Since r^A/d = 58.9, this choice satisfies our requirement, r^^ ^ L, and we take 
it as the fiducial size of the computational domain (see Table |l]). 

In these rescaled coordinates, Eqs. (Al) and (|A3[) are further transformed according to 



d(b ,d(j) d'^d) ,')d'^(b , ,,06 , . „, 

^7^. + (A6) 



dr dp' dr^ dp^ dp' 
dz d(' dz"^ (9(^2 Qf^' drdz dpd( 



= 5^. ?t = B'?i+BB'^4^ ?^-ABpt.. (A7) 



3. Discretization 



We discretize the computational domain as Pi — [i + l/2)h and Q — [j — N)h with i — 0, 1, . . . , iV and j 
0, 1, . . . , 2N where h is the spatial interval. Note that since 



L^{Nh) + ^{Nh)\ 

the fiducial parameter choices of Tab. [ijare achieved with a = 0.1 and A'' = 1600. 

The derivatives with respect to p and C are approximated by the central finite differences on the grid as 



(A8) 
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(j){pi,Q). The reason why we shift the p coordinate by half a spatial interval h is to 



be able to easily impose the Neumann boundary condition at p = 0. Considering an auxiliary point at i 
boundary condition reads 



dp 







0O,j - '/'-Ij 



0, 



90 J ■ 
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(AlO) 



p=0 



Thus the boundary condition for the value of the auxiliary point enforces the vanishing of the derivative without 
modifying the finite difference scheme. For example, the radial term in the Laplacian evaluated at z = becomes 
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Equation ( Al I written in (p, C) coordinate can then be recast as a matrix equation, schematically given as 

(linear terms of (/)a) = A^[0A, 0a] + A^[0A, 0b] — > = b = f (x) + b, (A12) 



where the matrix A represents the discrete linear operators on (j) given by the combination of Eq. ( A9 ) , being an 
asymmetric 9-band sparse (N + 1){2N + 1) x {N + 1){2N + 1) matrix, and b and f are vectors containing N[(j>A, 4>b] 
and [0a,0a] evaluated at (piXj)- 
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Iteration scheme 



To solve the set of nonlinear equations given in Eq. (A12), we use a combination of preconditioned Conjugate 



Gradient Squared (CGS) and Successive Over-Relaxatio n (SO R) methods according to Ref. [3T]. Let the n-th iterated 
solution be x*^"^. We evaluate the right-hand side of Eq. (A12|, b^") = b-f f(x(")). Then we solve the linear equation, 



Ax* = b("), by CGS. Finally, we update the solution x("+i) by SOR as 

x("+i) =cjx* + (1 -cj)x("\ 



(A13) 



where we set lu = 0.5 ~ 1 to make the iterative solution converge. 

To achieve the fast convergence of CGS, we precondition the matrix A with modified incomplete LU (MILU) 
decomposition, referred to as MILUCGS in the literature. We decompose A such that A — LDU + R, where L, U 
are just the copies of lower and upper triangular part of A, and D is a diagonal matrix, which are given by 






i < j 








i = j 


Aij 


i > j 




D, 



(AM) 



The diagonal element di is recursively calculated by 

d7^ = A 



Ai^i-ldi^l{Ai^l^i + j4i_i,i+Ar_i), 



Multiplying {LDU)~'^ in both sides in Eq. ( A12 1, we obtain a new operator matrix A' = {LDU)~'^A = I+{LDU)~'^R 
where we formally separated as A = LDU + R. The resultant matrix A' becomes close to a unit matrix in the sense 
that the weight of the non-zero components of A' becomes significant near the diagonal. In other words, the condition 
number k{A') = \\A'~'^\\/\\A'\\ becomes small. This treatment is frequently used to accelerate the convergence of 
CGS. 



We stop the SOR iteration in Eq. (A13I if we achieve 

||g(n+l) _ ylx("+^)| 



< esoR, 



where 



||b("+i)|| 

represents the standard 2-norm, and we set esOR = 10^® 



(A15) 



5. Convergence Tests 

First, we show the accuracy of our numerical results when we take the fiducial choice of parameters. In Fig. [8] we 



plot the ratio of the residuals of the field equation ( 26 ) to its source term 



2N[(j)A,4>B\ 

Setting esoR = 10"®, we find that the local residuals along z = —d are suppressed to less than 10"^ except in the 
far-field regime where they are still below lO"''. Recall that the boundary conditions at the edges of the box are set 
to enforce superposition. Given finite computational resources there is a trade-off between increased box size and 
central resolution. 

We therefore also test robustness of our results to the box size L and resolution h. In Fig. |9] (top), we vary a at 
fixed h OT N thus changing L. As long as L ^ r^^ results near the body are independent of box size. If L > r^j^ 
then we regain the far-field behavior, justifying the use of superposition boundary conditions. In Fig. [9] (bottom), we 
study resolution h by changing N at fixed a. Interior to body B, a resolution of h/rss < 0.125 is required to obtain 
converged results. 
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FIG. 8: The fractional accuracy of the numerical results with the fiducial choice of parameters along z = —d. 
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